In [1]:
# imports..
from sklearn.datasets import load_boston
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
import warnings
warnings.simplefilter('ignore')
In [2]:
# Gather Data....and Explore
boston_dataset=load_boston()
In [3]:
type(boston_dataset)
Out[3]:
In [4]:
boston_dataset
Out[4]:
In [5]:
dir(boston_dataset)
Out[5]:
In [6]:
print(boston_dataset.DESCR)
In [7]:
print(boston_dataset.data)
In [8]:
print(boston_dataset.feature_names)
In [9]:
print(boston_dataset.filename)
In [10]:
print(boston_dataset.target)
In [11]:
type(boston_dataset.data)
Out[11]:
In [12]:
boston_dataset.data.shape
Out[12]:
In [13]:
data = pd.DataFrame(data = boston_dataset.data, columns = boston_dataset.feature_names)
In [14]:
data
Out[14]:
In [15]:
data['PRICE'] = boston_dataset.target
In [16]:
data.head()
Out[16]:
In [17]:
data.tail()
Out[17]:
In [18]:
data.count()
Out[18]:
In [19]:
pd.isnull(data)
Out[19]:
In [20]:
pd.isnull(data).any()
Out[20]:
In [21]:
data.info()
In [22]:
plt.hist(data['PRICE'])
plt.show()
In [23]:
plt.hist(data['PRICE'], bins = 30)
plt.show()
In [24]:
plt.figure(figsize= (10,6))
plt.hist(data['PRICE'], bins =50, ec = 'black', color = '#800000', alpha = 0.7)
plt.xlabel('Price in 1000\'s')
plt.ylabel('Nr of Houses')
plt.show()
In [25]:
sns.distplot(data['PRICE'])
plt.show()
In [26]:
plt.figure(figsize= (10,6))
sns.distplot(data['PRICE'], color = 'red', bins = 50)
plt.show()
In [27]:
plt.figure(figsize= (10,6))
sns.distplot(data['PRICE'], color = 'red', bins = 50, hist=False)
plt.show()
In [28]:
plt.figure(figsize= (10,6))
sns.distplot(data['PRICE'], color = 'red', bins = 50,kde=False )
plt.show()
In [29]:
plt.figure(figsize= (10,6))
plt.hist(data['RM'], ec = 'black', color = 'pink')# removed bins
plt.xlabel('Average Number of Rooms')
plt.ylabel('Nr of Houses')
plt.show()
In [30]:
data['RM'].mean()
Out[30]:
In [31]:
plt.figure(figsize= (10,6))
plt.hist(data['RAD'], ec = 'black', color = 'green')
plt.xlabel('Accessibility to Highways')
plt.ylabel('Nr of Houses')
plt.show()
In [32]:
data['RAD'].value_counts()
Out[32]:
In [33]:
plt.figure(figsize= (10,6))
plt.hist(data['RAD'], ec = 'black', bins = 24, color = 'green')
plt.xlabel('Accessibility to Highways')
plt.ylabel('Nr of Houses')
plt.show()
In [34]:
accessibility = data['RAD'].value_counts()
In [35]:
accessibility.values
Out[35]:
In [36]:
print(accessibility.index)
In [37]:
type(accessibility)
Out[37]:
In [38]:
plt.figure(figsize= (10,6))
plt.bar(accessibility.index, accessibility )
plt.show()
In [39]:
plt.figure(figsize= (10,6))
plt.bar(accessibility.index, accessibility, color ='brown', ec = 'blue', width = 0.5)# changing the width, color and ec
plt.xlabel('Accessibility index to Highways')
plt.ylabel('Nr of Houses')
plt.show()
In [40]:
data['CHAS'].value_counts()
Out[40]:
In [41]:
data['PRICE'].min()
Out[41]:
In [42]:
data['PRICE'].max()
Out[42]:
In [43]:
data['PRICE'].mean()
Out[43]:
In [44]:
data['PRICE'].median()
Out[44]:
In [45]:
data.min()
Out[45]:
In [46]:
data.max()
Out[46]:
In [47]:
data.mean()
Out[47]:
In [48]:
data.median()
Out[48]:
In [49]:
data.describe()
Out[49]:
In [50]:
# Correlation
$$ \rho_{xy} = Corr(X,Y) $$¶
$$ -1.0 \leq \rho_{xy} \leq +1.0 $$¶
In [51]:
data['PRICE'].corr(data['RM'])
Out[51]:
In [52]:
data['PRICE'].corr(data['PTRATIO'])
Out[52]:
In [53]:
data.corr()
Out[53]:
In [54]:
# Heatmap
In [55]:
type(data.corr())
Out[55]:
In [56]:
plt.figure(figsize =(16,10))
sns.heatmap(data.corr(), annot = True, annot_kws ={'size':14} )
plt.show()
mask : boolean array or DataFrame, optional
If passed, data will not be shown in cells where mask is True.
Cells with missing values are automatically masked.
In [57]:
masker = np.zeros_like(data.corr())
In [58]:
masker
Out[58]:
In [59]:
type(masker)
Out[59]:
In [60]:
triangle_indices = np.triu_indices_from(masker)
In [61]:
triangle_indices
Out[61]:
In [62]:
masker[triangle_indices] = True
In [63]:
masker
Out[63]:
In [64]:
plt.figure(figsize =(16,10))
sns.heatmap(data.corr(), annot = True, annot_kws ={'size':14}, mask = masker )
plt.xticks(fontsize = 10)
plt.yticks(fontsize = 10)
plt.show()
In [65]:
nox_dis_corr = round(data['NOX'].corr(data['DIS']), 3)
plt.figure(figsize = (16, 8))
plt.title (f'Dist from Employment Centres vs Pollution , Corr ({nox_dis_corr})',
fontsize = 25, color = 'green')
plt.scatter(x=data['DIS'], y =data['NOX'], alpha = 0.5, color = 'red', s= 80)
plt.xlabel( 'Distance from Employment Centres...DIS', fontsize =14, color = 'green' )
plt.ylabel( 'Level of Pollutants Nitric Oxide...NOX', fontsize =14, color = 'green' )
plt.show()
In [66]:
sns.set()
sns.set_style('dark')
sns.jointplot(x=data['DIS'], y =data['NOX'], color = 'red', joint_kws = {'alpha' : 0.5},)
plt.show()
In [67]:
sns.set() # to reset any previous styling
sns.set_style('dark')# to set style.
plt.figure(figsize = (16, 8))
sns.scatterplot(x=data['DIS'], y =data['NOX'], color = 'red', alpha = 0.5) # joint plot creates scatter plot
plt.show()
In [68]:
sns.set()
sns.set_style('dark')
sns.jointplot(x=data['DIS'], y =data['NOX'], kind = 'hex' ) # kind is hex, darker hexagons in regions of higher density
plt.show()
In [69]:
sns.set() # to reset any previous styling
sns.set_style('dark')# to set style.
sns.jointplot(x=data['TAX'], y =data['RAD'], color = 'darkred', joint_kws = {'alpha' : 0.5} ) # jointplot gives scatter plot
plt.show()
In [70]:
sns.lmplot(x ='TAX', y = 'RAD', data = data)
plt.show()
In [71]:
sns.set() # to reset any previous styling
sns.set_style('dark')# to set style.
# height = 7 is optional
sns.jointplot(x=data['RM'], y =data['PRICE'], color = 'indigo', kind = 'scatter', height = 7) # jointplot gives scatter plot
plt.show()
In [72]:
sns.lmplot(x ='RM', y = 'PRICE', data = data, height =10, aspect = 2)
plt.show()
In [73]:
sns.lmplot(x ='RM', y = 'PRICE', data = data, height =10, aspect = 2, ci = None)
plt.show()
In [74]:
rm_price_corr = round(data['RM'].corr(data['PRICE']), 3)#calculating the corelation b/n NOX and DIS, also rounding to 3 decimals
plt.figure(figsize = (16, 8))
plt.title (f'Rooms vs Price , Corr ({rm_price_corr})',
fontsize = 25, color = 'green' )# using the f string notation to print the Correlatiin value in title
plt.scatter(x=data['RM'], y =data['PRICE'], alpha = 0.5, color = 'red', s= 80)
plt.xlabel( 'RM ....No of ROOMs ', fontsize =14, color = 'green' )
plt.ylabel( 'House price...PRICE', fontsize =14, color = 'green' )
plt.show()
In [75]:
# sns.pairplot(data)
# plt.show()
In [76]:
# sns.pairplot(data, kind = 'reg')
# plt.show()
Regression Multivariable¶
In [77]:
prices = data['PRICE'] # Target
features = data.drop('PRICE', axis = 1 ) # Features
In [78]:
X_train, X_test, y_train, y_test = train_test_split(features, prices, test_size = 0.2, random_state = 10)
In [79]:
len(X_train)/len(features)
Out[79]:
In [80]:
X_train.head()
Out[80]:
In [81]:
regr = LinearRegression()
regr.fit(X_train, y_train)
Out[81]:
In [82]:
np.round(regr.coef_, decimals =4)
Out[82]:
In [83]:
regr.intercept_
Out[83]:
In [84]:
pd.DataFrame(data = regr.coef_, index = X_train.columns, columns =['Coef'] )
Out[84]:
In [85]:
regr.score(X_train, y_train)
Out[85]:
In [86]:
regr.score(X_test, y_test)
Out[86]:
In [87]:
## Normal Distribution
In [88]:
from scipy.stats import norm
In [89]:
plt.figure(figsize = (10,5))
x_axis = np.arange(-10, 10, 0.001)
# mean = 0, SD = 3
plt.plot(x_axis, norm.pdf(x_axis,0,3))
plt.show()
In [90]:
# Data Transformations
In [91]:
plt.figure(figsize= (16,8))
plt.hist(data['PRICE'], bins =50, ec = 'black', color = 'blue', alpha = 0.5)# bins 50...setting color by hexcode
plt.xlabel('Price in 1000\'s')
plt.ylabel('Nr of Houses')
plt.show()
In [92]:
data['PRICE'].skew()
Out[92]:
In [93]:
y_log = np.log(data['PRICE'])
In [94]:
y_log.skew()
Out[94]:
In [95]:
sns.set()
plt.figure(figsize=(20,6))
sns.distplot(y_log)
plt.title(f'Log price: The skew is {round(y_log.skew(),2)}', fontsize = 25)# note the f string used to print the y_log.skew()
plt.xlabel('LOG PRICES ', fontsize =14, color = 'green' )
plt.show()
In [96]:
transformed_data = features # new data frame
transformed_data['LOG_PRICE'] = y_log # added the log price to the new data frame called transformed data
transformed_data
Out[96]:
In [97]:
sns.lmplot(x='LSTAT', y='PRICE', data=data, height=10, aspect =2,
scatter_kws={'alpha': 0.6}, line_kws={'color':'darkred'})
plt.show()
In [98]:
sns.lmplot(x='LSTAT', y='LOG_PRICE', data=transformed_data, height =10,aspect = 2,
scatter_kws={'alpha': 0.6}, line_kws={'color':'blue'})
plt.show()
In [99]:
prices = np.log(data['PRICE']) # Use log prices
features = data.drop('PRICE', axis=1)
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
regr = LinearRegression()
regr.fit(X_train, y_train)
print('Training data r-squared:', regr.score(X_train, y_train))
print('Test data r-squared:', regr.score(X_test, y_test))
pd.DataFrame(data=regr.coef_, index=X_train.columns, columns=['coef'])
Out[99]:
In [100]:
print('Intercept', regr.intercept_)
In [101]:
# Evaluating P values using statsmodel....and also the coefficients using the stats model
In [102]:
#------------------------------------------Comments-----------------------------------------------------------------------------
# 1. Using statsmodel.... please import ....import statsmodels.api as sm.....in the uppermost cell
# 2. Using statsmodel we can Evaluate, the p-values and coefficients
# 3. OLS stands for ORDINARY LEAST SQUARES, gives a linear regression model
# 4. We pass the y_train and X_train, into sm.OLS
# 5. We have to use sm.add_constant(X_train),only then we will be able to get the value of the constant(intercept)
# 6. We then create a pandas DataFrame, to hold the Coefficients and P-values
# 7. Deprecation warnings if any, maybe ignored or acted upon as in warnings, see import statement topmost cell
#------------------------------------------------------------------------------------------------------------------------------#
#results = sm.OLS(y_train, X_train.assign(const=1)).fit()
results = sm.OLS(y_train, sm.add_constant(X_train)).fit()# obtaining the regression coeffs including the constant(intercept).
print(results.params, '\n')# printing the values of the coefficients including the consstant(intercept)
print(results.pvalues,'\n')# printing the p values
pd.DataFrame({'coef': results.params, 'p-value': round(results.pvalues, 3)})# created DataFrame,p-values rounded to 3 decimals
# any p value more than 0.05 is not significant, look at the INDUS feature and AGE feature....
Out[102]:
In [103]:
# Multicollinearity
In [104]:
# 1. import.....from statsmodels.stats.outliers_influence import variance_inflation_factor
# 2. exog = using X_train with added constant.and getting the values from it.
# 3. sm.add_constant(x_train) will give a DataFrame
# 4. If you want to check do check (type(sm.add_constant(X_train)))
# 4. exog has to be a ndarray
# 5. Hence we used the .values, to passs only the values from the DataFrame.
# 3. exog.idx = index of the feature that we need to test
In [105]:
variance_inflation_factor(exog=sm.add_constant(X_train).values, exog_idx=10)# TAX is at index 10
Out[105]:
write a for loop that prints out all the VIFs for all the features¶
In [106]:
X_incl_const = sm.add_constant(X_train)
X_incl_const.shape[1]
Out[106]:
In [107]:
vif = []
for i in range(X_incl_const.shape[1]):
vif.append(variance_inflation_factor(exog=sm.add_constant(X_train).values, exog_idx=i))
print(vif)
In [108]:
pd.DataFrame({'coef_name': X_incl_const.columns, 'vif': np.around(vif, 2)})
Out[108]:
Bayesian Information Criterion .....the lesser the happier¶
Original model with log prices and all features ¶
In [109]:
X_incl_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_incl_const)
results_1 = model.fit()
org_coef = pd.DataFrame({'coef': results_1.params, 'p-value': round(results_1.pvalues, 3)})
print('BIC with all features and the log price is :' , results_1.bic, '\n')
print('r-squared is with all features and log price is :' , results_1.rsquared, '\n')
print(org_coef)
Model without INDUS¶
In [110]:
X_incl_const = sm.add_constant(X_train)
X_incl_const = X_incl_const.drop(['INDUS'], axis =1)
model = sm.OLS(y_train, X_incl_const)
results_2 = model.fit()
coef_minus_indus = pd.DataFrame({'coef': results_2.params, 'p-value': round(results_2.pvalues, 3)})
print('BIC without INDUS features and the log price is :' , results_2.bic, '\n')
print('r-squared without INDUS feature and log price is :' , results_2.rsquared, '\n')
print(coef_minus_indus)
Model without INDUS AND AGE¶
In [111]:
X_incl_const = sm.add_constant(X_train)
X_incl_const = X_incl_const.drop(['INDUS', 'AGE'], axis =1)
model = sm.OLS(y_train, X_incl_const)
results_3 = model.fit()
coef_minus_indus_age = pd.DataFrame({'coef': results_3.params, 'p-value': round(results_3.pvalues, 3)})
print('BIC without INDUS features and the log price is :' , results_3.bic, '\n')
print('r-squared without INDUS feature and log price is :' , results_3.rsquared, '\n')
print(coef_minus_indus_age)
In [112]:
X_incl_const
Out[112]:
In [113]:
X_incl_const = sm.add_constant(X_train)
In [114]:
X_incl_const
Out[114]:
In [115]:
X_incl_const.shape
Out[115]:
In [116]:
X_train.shape
Out[116]:
RESIDIAL PLOT OF MODEL 3 THE LOG MODEL MINUS INDUS AND AGE¶
In [117]:
plt.figure(figsize=(20,6))
plt.scatter(x=results_3.fittedvalues, y=results_3.resid, c='brown', alpha=0.8)
plt.xlabel('Predicted log prices $\hat y _i$', fontsize=20)
plt.ylabel('Residuals', fontsize=20)
plt.title('Residuals vs Fitted Values of the reduced log model', fontsize=20)
plt.show()
# Mean Squared Error & R-Squared
reduced_log_mse_minus_indus_minus_age = round(results_3.mse_resid, 3)# for later use and comparision
reduced_log_rsquared_minus_indus_minus_age = round(results_3.rsquared, 3) # for later use and comparision
CREATING A CRAZY ...OMITTED VARIABLES MODEL: INDUS , AGE , LSTAT , RM, NOX, CRIM¶
In [118]:
X_incl_const = sm.add_constant(X_train)
X_incl_const = X_incl_const.drop(['INDUS', 'AGE','LSTAT', 'RM', 'NOX', 'CRIM' ], axis =1)
model = sm.OLS(y_train, X_incl_const)
results_4 = model.fit()
coef_minus_many_variables = pd.DataFrame({'coef': results_4.params, 'p-value': round(results_4.pvalues, 3)})
print('BIC after dropping many variables and the log price is :' , results_4.bic, '\n')
print('r-squared after dropping many variables and log price is :' , results_4.rsquared, '\n')
print(coef_minus_many_variables)
RESIDUAL PLOT: OMITTED VARIABLES MODEL 'INDUS', 'AGE','LSTAT', 'RM', 'NOX', 'CRIM'¶
In [119]:
plt.figure(figsize=(20,6))
plt.scatter(x=results_4.fittedvalues, y=results_4.resid, c='brown', alpha=0.8)
plt.xlabel('Predicted log prices after omitting many variables $\hat y _i$', fontsize=20)
plt.ylabel('Residuals', fontsize=20)
plt.title('Residuals vs Fitted Values after ommitting many variables', fontsize=20)
plt.show()
# Mean Squared Error & R-Squared
ommited_var_mse = round(results_4.mse_resid, 3)# for later use and comparision
ommitted_var_rsquared = round(results_4.rsquared, 3) # for later use and comparision
GOING BACK TO THE REDUCED LOG MODEL WITHOUT INDUS AND AGE¶
In [120]:
resid_mean = round(results_3.resid.mean(), 3)
resid_mean
Out[120]:
In [121]:
resid_skew = round(results_3.resid.skew(), 3)
resid_skew
Out[121]:
In [122]:
plt.figure(figsize= (20,6))
sns.distplot(results_3.resid, color='navy')
plt.title(f'Residual Dist of Log price model w/o INDUS and AGE:residuals Skew ({resid_skew}) Mean ({resid_mean})',fontsize = 20)
plt.show()
In [123]:
frames = [org_coef, coef_minus_indus, coef_minus_indus_age, coef_minus_many_variables]
pd.concat(frames, axis=1, sort = False)
Out[123]:
In [124]:
print(f'BIC WITH ALL:{round(results_1.bic)}:W/O INDUS {round(results_2.bic)}:W/O INDUS AND AGE {round(results_3.bic)}:W/O MANY VARIABLES {round(results_4.bic)}')
In [125]:
print(f'RSQD:WITH ALL {round((results_1.rsquared),2)}: W/O INDUS {round((results_2.rsquared),2)}: W/O INDUS,AGE {round((results_3.rsquared),2)}:W/O MANY VARIABLES {round((results_4.rsquared),2)}')
In [126]:
# Our Estimate for a house is 30000....calculate the upper bound and the lower bound using the log model_minus_indus_minus_age
sd_1 = round(np.sqrt(results_3.mse_resid),3)
sd_2 = 2*(round(np.sqrt(results_3.mse_resid),3))
print('1 s.d. in log prices is', sd_1)
print('2 s.d. in log prices is', sd_2)
upper_bound = np.log(30) + sd_2
print('The upper bound in log prices for a 95% prediction interval is ', upper_bound)
print('The upper bound in normal prices is $', np.e**upper_bound * 1000)
lower_bound = np.log(30) - sd_2
print('The lower bound in log prices for a 95% prediction interval is ', lower_bound)
print('The lower bound in normal prices is $', np.e**lower_bound * 1000)
In [127]:
data
Out[127]:
EXAMINING THE DOLLAR PRICE MODEL WITHOUT THE LOG CONVERSION¶
In [148]:
house_prices = data['PRICE']
In [149]:
house_features = data.drop(['PRICE'], axis =1)
In [150]:
X_train, X_test, y_train, y_test = train_test_split(house_features, house_prices,
test_size=0.2, random_state=10)
In [151]:
X_incl_const_to_house_features = sm.add_constant(X_train)
In [152]:
model = sm.OLS(y_train, X_incl_const_to_house_features)
In [153]:
results_5 = model.fit()
In [154]:
round(results_5.params,3)
Out[154]:
In [155]:
results_5.bic
Out[155]:
In [156]:
results_5.rsquared
Out[156]:
In [157]:
round(results_5.pvalues,2)
Out[157]:
In [158]:
plt.figure(figsize=(20,6))
plt.scatter(x=results_5.fittedvalues, y=results_5.resid, c='brown', alpha=0.8)
plt.xlabel('Predicted Original prices $\hat y _i$', fontsize=20)
plt.ylabel('Residuals', fontsize=20)
plt.title('Residuals vs Fitted Values of the original model', fontsize=20)
plt.show()
# Mean Squared Error & R-Squared
original_mse = round(results_5.mse_resid, 3)# for later use and comparision
original_rsquared = round(results_5.rsquared, 3) # for later use and comparision
In [159]:
plt.figure(figsize= (20,6))
sns.distplot(results_5.resid, color='navy')
plt.title(f'Residual Dist Dollar price model)',fontsize = 20)
plt.show()
In [160]:
print(f'Skew:({round(results_5.resid.skew(),2)}), residuals mean : ({round(results_5.resid.mean(),2)}))')
In [ ]:
In [ ]:
In [ ]: